#fit tadpoles model
tadpoles <- read.table("ecol 563/tadpoles.csv", header=T, sep=',')
tadpoles$fac3 <- factor(tadpoles$fac3)
out1 <- lm(response~fac1 + fac2 + fac3 + fac1:fac2 + fac1:fac3 + fac2:fac3 + fac1:fac2:fac3, data=tadpoles)
out2 <- lm(response~fac1 + fac2 + fac3 + fac1:fac2  + fac2:fac3, data=tadpoles)
anova(out2)
summary(out2)
 
# components of summary object
names(summary(out2))
# summary table
summary(out2)$coefficients
 
# extract effect estimates and std errs (omit intercept)
ests <- summary(out2)$coefficients[2:8,1]
std.errs <- summary(out2)$coefficients[2:8,2]
 
# degrees of freedom for confidence intervals
names(out2)
out2$df.residual
 
#confidence intervals
low95 <- ests+qt(.025, out2$df.residual)*std.errs
up95 <- ests+qt(.975, out2$df.residual)*std.errs
low95
up95
 
# can also obtain confidence intervals with confint function
confint(out2)
 
#assemble components in a single data frame
# make a factor in which main effects are listed first
new.data <- data.frame(var.labels=factor(names(ests), levels=names(ests)), ests, low95, up95)
new.data
levels(new.data$var.labels)
 
# in base graphics dotchart; in lattice graphics dotplot
library(lattice)
dotchart(new.data$ests, labels=new.data$var.labels, xlab='Effect estimates')
dotplot(var.labels~ests, data=new.data, xlab='Effect estimates')
 
# Effects graph using lattice graphics
 
dotplot(var.labels~ests, data=new.data, xlab='Effect estimates', xlim=c(min(new.data$low95)-.02, max(new.data$up95)+.02), panel=function(x,y){
# use dotplot to get gridlines
panel.dotplot(x,y,col='white')
# add confidence intervals
panel.segments(new.data$low95, as.numeric(y), new.data$up95, as.numeric(y))
# add points
panel.xyplot(x,y,pch=16,cex=1.2)
# add vertical line at zero
panel.abline(v=0, col=2, lty=2)
})
 
# Use same color for intervals and point estimates. Change bar ends to butted
dotplot(var.labels~ests, data=new.data, xlab='Effect estimates', xlim=c(min(new.data$low95)-.02, max(new.data$up95)+.02), panel=function(x,y){
 panel.dotplot(x, y, col='white')
 panel.segments(new.data$low95, as.numeric(y), new.data$up95, as.numeric(y), col="#0080ff", lineend=1)
 panel.xyplot(x,y,pch=16,cex=1.2)
 panel.abline(v=0, col=2, lty=2)
 })
 
# math expressions for y-axis labels
mylabel<-c('Normal', 'Ru486', 'Shrimp', 'Sibship', expression('Normal'%*%'Shrimp'), expression('Ru486'%*%'Shrimp'), expression('Shrimp'%*%'Sibship'))
mylabel
 
# change labels on y-axis
dotplot(var.labels~ests, data=new.data, xlab='Effect estimates', xlim=c(min(new.data$low95)-.02, max(new.data$up95)+.02), panel=function(x,y){
 panel.dotplot(x, y, col='white')
 panel.segments(new.data$low95, as.numeric(y), new.data$up95, as.numeric(y), col="#0080ff", lineend=1)
 panel.xyplot(x, y, pch=16, cex=1.2)
 panel.abline(v=0, col=2, lty=2)
 }, scales=list(y=list(at=1:7, labels=mylabel)))
 
# demonstration of mathematical typesetting in R
demo(plotmath)
 
# Effects graph using base graphics
 
# x-axis limits
range(new.data[,2:4])
# set limits for x-axis and y-axis
plot(range(new.data[,2:4]), c(1-.3,7+.3))
#set up graph but without plotting
plot(range(new.data[,2:4]), c(1-.3,7+.3), type='n', xlab='Effect estimates', ylab='',axes=F)
# x-axis
axis(1)
# y-axis
axis(2, at=1:7, labels=mylabel, las=2)
box()
 
# increase margins on left
par("mar")
par("mar")->oldmar
par(mar=c(5.1,8.1,4.1,2.1))
 
# redo graph
plot(range(new.data[,2:4]), c(1-.3,7+.3), type='n', xlab='Effect estimates', ylab='',axes=F)
axis(1)
axis(2, at=1:7, labels=mylabel, las=2)
box()
 
# create a numeric variable for the labels
new.data$num.labels <- as.numeric(new.data$var.labels)
new.data
# change bar ends to butted (from default round)
par(lend=1)
# confidence intervals
segments(new.data$low95, new.data$num.labels, new.data$up95, new.data$num.labels, col="#0080ff")
# point estimates
points(new.data$ests, new.data$num.labels, pch=16, cex=1.2, col="#0080ff")
# vertical line at zero
abline(v=0, lty=2, col=2)
# add grid lines to graph
grid(nx=NA, ny=NULL, col="grey70")
 
#reset margins
par(mar=oldmar)
 
#### Obtaining treament means and their confidence intervals #####
 
## Method 1: use the effects package
 
library(effects)
effect('fac1:fac2', out2)
summary(out2)
ef1a <- effect('fac1:fac2', out2, given.values=c(fac32=0), se=T)
ef2a <- effect('fac1:fac2', out2, given.values=c(fac32=1), se=T)
names(ef1a)
 
#treatment means
ef1a$fit
# standard errors of means
ef1a$se
 
# 95% confidence intervals
ef1a$lower
ef1a$upper
 
plot(ef1a,multiline=T)
 
names(ef1a)
# factor levels
ef1a$x
 
# mean statistics for sibship 1
part1 <- data.frame(ef1a$x, fac3=1, est=ef1a$fit, se=ef1a$se,low95=ef1a$lower, up95=ef1a$upper)
part1
# mean statistics for sibship 2
part2 <- data.frame(ef2a$x, fac3=2, est=ef2a$fit, se=ef2a$se,low95=ef2a$lower, up95=ef2a$upper)
part2
# combine both sibships in a single matrix
part.eff <- rbind(part1, part2)
part.eff
 
# Method 2: obtain means and std errs using predict function
fac.vals2 <- expand.grid(fac1=c('Co','No','Ru'), fac2=c('D','S'), fac3=1:2)
fac.vals2
 
#we need to convert fac3 to a factor or we get an error message
out.p <- predict(out2, newdata=fac.vals2, se.fit=T, interval="confidence")
fac.vals2$fac3 <- factor(fac.vals2$fac3)
out.p <- predict(out2, newdata=fac.vals2, se.fit=T, interval="confidence")
names(out.p)
#estimated means and CIs
out.p$fit
#standard error of means
out.p$se.fit
results.pred <- data.frame(fac.vals2, out.p$fit, out.p$se.fit)
results.pred
names(results.pred)[4:7] <- c('est', 'low95', 'up95', 'se')
results.pred
coef(out2)
 
 
# Method 3: obtain means and std errs by hand
 
#function to create dummy vector 
myvec <- function(fac1, fac2, fac3) c(1, fac1=='No', fac1=='Ru', fac2=='S', fac3==2, (fac1=='No')*(fac2=='S'), (fac1=='Ru')*(fac2=='S'), (fac2=='S')*(fac3==2))
myvec('Co','D',1)
myvec('Co','S',1)
# means for different treatments
myvec('Co','D',1) %*% coef(out2)
myvec('No','S',1) %*% coef(out2)
fac.vals2
#obtain dummy vectors for each treatment
apply(fac.vals2, 1, function(x) myvec(x[1],x[2],x[3]))
#make dummy vectors the rows of a matrix
t(apply(fac.vals2, 1, function(x) myvec(x[1],x[2],x[3])))
 
#contrast matrix
cmat <- t(apply(fac.vals2, 1, function(x) myvec(x[1],x[2],x[3])))
 
#treatment means
ests <- cmat%*%coef(out2)
ests
 
#variance-covariance matrix of treatment means
vmat <- cmat %*% vcov(out2) %*% t(cmat)
dim(vmat)
 
# std errs of means
se <- sqrt(diag(vmat))
se
 
# 95% confidence intervals
low95 <- ests + qt(.025,out2$df.residual)*se
up95 <- ests + qt(.975,out2$df.residual)*se
 
#assemble results in data frame
results.mine <- data.frame(fac.vals2, ests, se, low95, up95)
results.mine
 
#create a numeric version of the fac1 for plotting
results.mine$fac1num <- as.numeric(results.mine$fac1)
 
#restrict data to sibship 1
results.mine.sub <- results.mine[results.mine$fac3==1,]
results.mine.sub
 
#set up axes for interaction graph
plot(ests~fac1num, data=results.mine.sub, xlim=c(0.8,3.2), ylim=range(results.mine.sub[,6:7]), ylab='Mean mitotic activity', xlab='Hormonal treatment', type='n', axes=F)
axis(2)
axis(1, at=1:3, labels=levels(results.mine.sub$fac1))
box()
 
# shift amount for plotted profiles
jitter <- 0.05
 
#Detritus profile
cur.dat<-results.mine.sub[results.mine.sub$fac2=='D',]
points(cur.dat$fac1num-jitter, cur.dat$ests, col=1, pch=16)
lines(cur.dat$fac1num-jitter, cur.dat$ests, col=1)
arrows(cur.dat$fac1num-jitter, cur.dat$low95, cur.dat$fac1num-jitter, cur.dat$up95, col=1, code=3, angle=90, length=.05)
 
#Shrimp profile
cur.dat<-results.mine.sub[results.mine.sub$fac2=='S',]
points(cur.dat$fac1num+jitter, cur.dat$ests, col=2, pch=15)
lines(cur.dat$fac1num+jitter, cur.dat$ests,col=2, lty=2)
arrows(cur.dat$fac1num+jitter, cur.dat$low95, cur.dat$fac1num+jitter, cur.dat$up95, col=2, code=3, angle=90, length=.05)
 
#legend and margin text
legend('topleft', c('Detritus','Shrimp'), col=c(1,2), pch=c(16,15), lty=c(1,2), pt.cex=1, cex=.9, bty='n', title='Food type')
mtext('Sibship 1', side=3, line=.5, cex=.9)